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ABSTRACT 


This thesis deals with a simulation of a missile versus target engagement 
scenario. After deriving simplified transfer functions for the missile seeker 
head, missile autopilot, missile dynamics, and target dynamics, a three 
dimensional simulation is developed using classical proportional navigation. The 
scenario is simulated using state variable design. A forward time solution of the 
two dimensional problem is developed which is converted to an adjoint model. 
The adjoint model is used to determine the optimal time to initiate simplified and 
tactical evasion maneuvers in order to maximize the final miss distance. 








THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this 
research may not have been exercised for all cases of interest. While every 
effort has been made, within the time available, to ensure that the programs are 
free of computational and logical errors, they cannot be considered validated. 
Any application of these programs without additional verification is at the risk of 
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I. INTRODUCTION 


Since the early years of the Vietnam War the surface-to-air missile has 
become the largest threat to combat aviators. During the emergent age of low 
intensity conflicts, a goal of zero percent aircraft losses must remain a high 
priority to strike planners and decision makers. In order to minimize losses to 
enemy antiair weapon systems, in particular surface-to-air missiles, these 
systems must be generally understood and tactics necessary to defeat these 
missiles must be investigated. 

Aircraft combat survivability is defined as "the capability of an aircraft to 
avoid and/or withstand a man-made hostile environment." [Ref. 1] Figure 1 
displays the basic aircraft survivability relationships. The probability of an 
aircraft being hit by antiair weapons is termed Ph which is referred to as the 
susceptibility of the aircraft. Susceptibility is calculated by the product 

p h~ p a p Drr p LGD (i* 1 ) 

where i 

Pa = probability the threat is active and ready to engage the 

aircraft 

Pdtt * probability the aircraft is detected, identified, and tracked. 

Plod = probability the threat is launched, guided and either hits the 

aircraft or detonatesdose enough to cause a hit 

The vulnerability of an aircraft is defined as the inability to withstand the 

damage caused by a hostile environment It can be measured by Pjc/h which is 

the probability of an aircraft kill given it is hit by hostile fire. This leads to the 

relationship 

probability of kill - susceptibility-vulnerability, (1.2) 

1 







Figure 1. Aircraft Survivability Diagram 


or in other words, 

Pk-Pk'Pic/H. 0*3) 

It is obvious then that the probability of survival is 

P S »1-P K . [Ref. 1] (1.4) 

Therefore. in order to increase the survivability of an aircraft, the probability 
of kill must be reduced. One way that this can be achieved is by decreasing 
Plcd* The probability that a missile is successfully guided and either hits an 
aircraft or detonates close enough to hit it, can be reduced by degrading the 
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guidance and control system of the missile. The employment of high 
acceleration turns and random maneuvers are examples of methods that can be 
used. 

This research develops a missile/target simulation program using 
classical proportional navigation. A three dimensional model is produced 
assuming a dual gimbal axis seeker head. Chapter II introduces the idea of 
proportional navigation and the actual guidance law is developed. In Chapter 
HI the transfer functions for the individual missile subsystems are determined. 
The specifics of the computer simulation geometry are discussed leading to the 
actual relationships used in the program. Additionally, a two dimensional 
forward time model is developed which is converted to an adjoint model. The 
adjoint model is a useful tool for analyzing time varying systems. It is used in 
this thesis to determine optimal miss distance parameters. Chapter W consists 
of actual simulation results. Two typical scenarios are conducted to show the 
capabilities of the program against various target maneuvers. The adjoint 
model is then used to determine an optimal time to initiate various evasion 
maneuvers. Conclusions and recommendations follow in Chapter V. 

All computer simulations are developed and conducted using the Matrix 
Laboratory (MATLAB) and the three dimensional plots are generated using 
die Display Integrated Software System and Plotting Language (DISSPLA). 
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II. BASIC PROPORTIONAL NAVIGATION 


A. GENERAL 

Proportional navigation is the basic guidance law used in the majority of the 
threat missiles in operational use today. This method of guidance generates 
missile acceleration commands proportional to the line of sight (LOS) rate. 
Figure 2 shows the basic parameters and geometry associated with the two 
dimensional missile/target engagement problem and the parameters are defined 
below: 

V| = target velocity. 

V m = missile velocity. 

7t = target flight path angle. 

Ym = missile flight path angle. 

R = missile to target range. 

Rt = target range. 

R m = missile range. 

o = line of sight angle. 

Ct = target line of sight angle. 

Om = missile line of sight angle. 

B. CONSTANT BEARING COURSE 

A constant bearing course is one where the line of sight between the 
target and the missile maintains a constant orientation in space. As a result, 
progressive lines of sight remain parallel to each other as the engagement 
procedes. Similarly, the line of sight angle, o, remains constant. Figure 3 






ftafertnc* 


Figure 2. Missile/Target Geometry 

depicts the constant bearing course idea. As long as there is a positive closing 
velocity between the missile and the target, the constant bearing course 
concept will ensure an intercept. Proportional navigation uses the constant 
bearing course idea by driving the line of sight rate, a, to zero. [Ref. 2] 










Figure 3. Constant Bearing Course 


C. PROPORTIONAL NAVIGATION SCHEME 

Figure 4 depicts the basic proportional navigation scheme. Assuming that 
the seeker head of the missile follows the target, the transverse acceleration 
perpendicular to the lin , of sight will equal the acceleration of the R vector in 
that direction. Mathematically, the acceleration of R is 

A* = (R + <*> x to x R)i* + (2to x R x R+© x r)? # (2.1) 

where 

R = missile/target line of sight vector. 
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R = closing rate along R. 

R = acceleration along R. 

to = angular rate of change of R in inertial 

space. 

Am - missile acceleration perpendicular to 

A t = target acceleration perpendicular to R. 

Ar sb overall acceleration of R. 


At this point, a missile acceleration, A m , equal to the target acceleration. At, will 
make the line of sight parallel to its original direction. As long as R remains 



Figure 4, Vectorial Proportional Navigation Scheme 
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along R (<o=0) a missile/target impact is assured. So, the transverse 
acceleration command is 

A t -A m =<bxR+2(©xR). (2.2) 

Assuming the line of sight rate is equal to the angular rate of change of R in 
inertial space, equation (2.2) now becomes 

A t -A m =Ra+2Ra. (2.3) 

D. PROPORTIONAL NAVIGATION GUIDANCE LAW 

In the classical proportional navigation scheme, the missile course is one in 
which the rate of change of the missile heading is directly proportional to the 
rate of rotation of the line of sight vector from the missile to the target. As a 
result, this course change is intended to counteract the rotation of the line of 
sight, thus returning to a constant bearing course. The movement of the missile 
and target cause the line of sight to rotate resulting in a differential displacement 
between the missile and the targe, perpendicular to the range line. Figure 5 
depicts this geometry. [Ref. 3] The proportional navigation guidance law 
attempts to generate an acceleration command, A«, perpendicular to the line of 
sight. 

Assume a gyro stabilized seeker head. If there is no torque applied to the 
gyro, the seeker will not rotate. Assuming the seeker tracks; the target, the 
gyro angle will follow the line of sight Applying die equation of motion for a 
gyro stabilized seeker 

L*I<pQ (2.4) 

where 

L - applied torque, 

oa = Spin angular velocity. 
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I = moment of inertia of the gyroscope. 

Q = rate of precession of the gyroscope. 

Applying this to the case when the seeker head tracks the target, Q is then 
replaced by the rate at which the gyro is torqued in space. This is simply 6 
which is the line of sight rate. [Ref. 4] Thus, equation (2.4) becomes 

L = Ioxj. (2.5) 


This torque is in turn applied to the control surface of the missile leading to 
the relationship 

A m = kL = klcoa (2.6) 

where k is a constant of proportionality. Referring to Figure 6, a relationship is 
determined for Am in terms of the rate of change of the missile flight path angle, 
y m . Given the missile velocity vector at some point in time, V n (t), and suppose 
the missile undergoes an acceleration. Am. during an interval of time, dt The 
velocity vector is then displaced and is represented by the vector V m (t+dt). 
The angle the vector is traversed is simply dy m , the differential missile flight 
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path angle. [Ref. 4] For small angles (which are guaranteed by making dt 
small) the following relationship is obvious 

Amdt = Vmdym. (2.7) 

Dividing equation (2.7) by the time interval, dt, the missile acceleration is 

defined as 

Am * V m —jo. s V m y m . (2.8) 

Comtv „ Rations (2.6) and (2.8) 

V m Y ni *kI(ocr. (2.9) 

Dividing through by V m , the proportional navigation guidance law becomes 




( 2 . 10 ) 
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Equation (2.11) represents the classical proportional navigation equation where 

Y m = rate of change of the missile heading, 

o = rate of change of the line of sight. 

N = proportional navigation ratio. 

The navigation ratio determines the sensitivity of the missile system. A high 
navigation ratio will lead to rather high gains resulting in large missile 
commands for small changes in the line of sight rate. On the other hand, small 
values for N will lead to small missile commands for a given d. Larger 
navigation ratios are preferred for head on engagements and smaller ones are 
preferred for tail chase cases. For this research the navigation ratio between 
three and five was chosen [Ref. 4]. 





HI. SIMULATION DEVELOPMENT 


A. INTRODUCTION 

The development of the actual seeker head and autopilot transfer functions 
is discussed. Ideal missile and target trajectory equations are used and the 
state equations for these systems are generated. Basic geometric relationships 
are used to develop the missile and target proportional navigation equations. 
Once these equations are determined for the two dimensional problem, the state 
equations are then augmented to a three dimensional problem. The continuous 
state equations for the seeker head, autopilot, missile kinematics, and target 
kinematics are converted to equations appropriate for digital simulation. 
Finally, the two dimensional forward time model is developed and converted to 
an adjoint model. 

B. MISSILE SUBSYSTEM DEVELOPMENT 

1. Subsystem Layout 

A basic functional block diagram of a generic tactical missile system is 
shown in Figure 7. Although the exact configuration and description of each 
component depends on many different factors, an attempt will be made to 
develop a basic system sufficient for simplified simulation. Figure 8 shows the 
typical physical location of the various missile subsystems. 

2. Seeker Head Development 

The seeker head can simply be thought of as the eye of the missile. It 
is able to detect, acquire, and track a target by sensing some unique 
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Figure 7. Missile Subsystem Block Diagram 


characteristic of the target itself. This usually consists of the radiation or 
reflection of energy by the target. 

A seeker with a narrow field of view will be used. Figure 9 shows a 
basic gimballed seeker head configuration. Here, the actual seeker is mounted 
on a gim balled platform and it maintains the target within the field of view by 
rotating the platform. The inertial rotation rate of the line of sight provides the 
missile with the required tracking information. [Ref. 5] 

Figure 10 displays a diagram for an actual seeker head where 
|3 = seeker head gimbal angle. 

The control torque to the seeker results in the following equation of motion 

T*Ip (3.1) 

where 

T = applied torque. 

I * moment of inertia. 

3 * angular acceleration. 
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Figure 8. Basic Missile Layout 

Solving for p 

P-y~ ki(P-o)-k2p (3.2) 

a-k2p-kjp+k,^. (3.3) 

The constants ki and k 2 are determined by the time constant of the seeker head 
system. 





Now taking the Laplace Transform of equation (3.3) and assuming 


zero initial conditions 

£{P} = JC,{-kjMiP+M?} 


IS 


(3.4) 







s 2 P(s)»-k2sP(s)-kiP(s)+k 1 a(s) (3.5) 

8 2 p(s)+k 2 sp(s)+k 1 p(s) = k|a(s). (3.6) 

Finally, the transfer functior becomes 

. k l. .*1—(3.7) 

o(s) r + k 2 s+kj (s+l/t*) 2 

For this simulation the time constant of the seeker head, t*, was 
chosen to be 0.1 second which is a good approximation of a real world system. 
From equation (3.7) the constants ki and k 2 are defined in terms of this time 
constant where 

ki*(l/**) 2 -100, (3.8) 

k 2 *2(l/t*)«20. (3.9) 
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Figure 11 depicts the signal flow graph of the seeker head system. 
From this diagram the continuous time state equation of the form 

*sh = AXsh+BUjh (3.10) 

is easily determined. Selecting the state vector to be 


■CM 


and the input u^ = a, equation (3.10) becomes 


**“{-100 - 2 o] Xsh + [i Ush * 


(3.11) 


(3.12) 



Figure 11. Seeker Head Signal Flow Graph 
3. Guidance and Autopilot Development 

The guidance law within the missile system determines the best 
trajectory for the missile based on the missile position, target position, missile 
capability and the desired objectives. A command is sent to the autopilot which 
determines the control (i.e., actuator position and thrust) necessary to perform 
such a command. In the proportional navigation guidance law, the missile 
command* are generated in order to change the missile flight path rate 
proportionally to the missile to target line of sight rate. 
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In this simulation a very simplified autopilot will be used. Referring to 
Figure 12, the transfer function is easily developed. The angle of attack of the 
missile is assumed to be zero, so the velocity of the missile will be aligned with 
the missile center line at an angle, y m . Applying a torque to the missile about 
the center of gravity results in tl following equation of motion 

Tcon, = I C gY (3.13) 

where 

= commanded torque. 

leg = moment of inertia about the 

center of gravity. 

y = angular acceleration of the 

missile flight path. 



Figure 12. Basic Missile Layout 
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Solving for y. 


y=^ ffi -=-lcy n +kNp (3.14) 

leg 

where k is determined by the slowest time constant of the missile/autopilot. 

Now taking the Laplace Transform of equation (3.14) and assuming 
zero initial conditions 

X,{Ym}=/-{fcrm+kNp} 

S 2 Ym 00 * -fcYm (0 + kNsTn, (s) 

(s 2 + ks)y m (s) = kNsf5(s). 

Finally, the transfer function becomes 

Yn, (s) _ kN 
p(s) “s+k 

where k is defined above. The autopilot time constant ,t^,,was chosen to be 1.0 
second, so 

k = l/Tjp *1.0. (3.19) 

Figure 13 shows the signal flow graph for the missile autopilot system. 

From this diagram the continuous state equation is easily determined. Selecting 
the state vector to be 

VM-ftm] (3-20) 

and the input u^pNp, the state equation becomes 

at [“* 1 ] x V + [ 1 ] u ap- (3.21) 


(3.15) 

(3.16) 

(3.17) 

(3.18) 
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4. Missile and Target Dynamics 


The signal flow graphs for the target and missile dynamics are shown 
in Figure 14 where A m = V m y m . A t is selected independently. The state 
vectors are 


x m = 




y m 

Lym 


and 


(3.22) 


*t = 


y t 

Ly*J 


With u m * "j and u t = I , the state equations become 

L»myJ L u tyJ 


(323) 


20 







• 

Jn_ 

• • • 

v m 1/s 

1/s 

Vm 


Missile Dynamics 



A, 

i Vt vs h 

1/s 

y t 

Target Dynamics 


Figure 14. Missile and Target Dynamics 


and 


'0 10 0] TO O' 

0 0 0 0 1 0 

x ® = 0 0 0 1 X ® + 0 o u ® 

.0 0 o oj [o 1. 


(3.24) 


0 10 0 ' 
0 0 0 0 
0 0 0 1 
0 0 0 0 


(3.25) 


5. Overall System Layout 

Figure 15 shows the total system signal flow graph. The seeker head 
angle rate, {3, is the estimate of the line of sight rate, 6. 

6. Three Dimensional State Definitions 

Assuming that the missile is roll stabilized and controllable in both pitch 
and yaw, the state equations for the missile subsystems are extended. 
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Auto Pilot 


Missile Dynamics 








The seeker head state matrix now becomes 



Py*w 


(3.26) 


and the continuous state equation is obtained by simply augmenting the two 
dimensional problem. 

The resulting equation becomes 


Xjj, = 


' 0 

1 

0 

0 ' 

- 

* 0 O' 

-100 

-20 

0 

0 


100 0 

0 

0 

0 

1 1 

*sh + 

0 0 

0 

0 

-100 

-20 J 


. 0 10 °. 


u sh 


c\ 


where u^, 

The autopilot state matrix now becomes 


[ Ym_pitch 1 

t«.^w J 


(327) 


(328) 


and the continuous state equation is again simply obtained from an 
augmentation of the two dimensional equation. The resulting equation for the 
autopilot is 


=[o -,Mo l]"** 


(329) 
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where 


:i 


P pitch 
L Pyaw 


For the missile and target dynamics equations the state equations are 
simply augmented with vertical position and vertical velocity states. 

These equations are 


and the state variable dynamics become 


(330) 


'0 

1 

0 

0 

0 

or 


0 

0 

O' 

0 

0 

0 

0 

0 

0 


1 

0 

0 

0 

0 

0 

1 

0 

0 

x+ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 


0 

0 

0 

0 

0 

0 

0 

0 

0, 


0 

0 

1. 


(331) 


where u« 


L u * 


, the force per unit mass acting on the missile or target. 
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C. MISSILE AND TARGET GEOMETRY 
1. Velocity Relationships 

Figure 16 shows the vectorial relationship of the missile and target 
velocity. From this figure 

= ( 3 - 32 ) 

and 

V t =V V * 2 +y+V tt 2 . (333) 



In order to simulate the scenario in three dimensional space a line of 
sight angle in pitch and yaw must be defined. Figure 17 depicts these angular 
relationships. Given the target and missile position in Cartesian coordinates the 
line of sight angles are defined as 

<*pfch * tan ~ l [(*t -*■)/'^( x t“*m ) 2 + (y» - 71 f ] (3-34) 
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Figure 17. Pitch and Yaw Line of Sight Angle Definitions 


and 


O y . w = tan" , [(y t -y m )/(x t -x ni )J. (335) 

Figure 18 shows the target and missile line of sight ang>e definitions and these 
angles are expressed mathematically as 

°m_pad> = tan -1 [z m /jx a 2 + y m 2 ], (336) 

^n_ytw =tan“ , [y m /x ni ], (337) 

and 


a t_piich * tan” 1 [*t / V x t 2+ yt 2 ]. (338) 

<*»_*•** tan -1 [y,/x t ]. ( 339 ) 





26 






z 


malortgt (t.y.x) 



Figure 18. Target and Missile Line of Sight Angle Definitions 
3. Flight Path Angles 


The missile and target flight path angles are also defined in three 
dimensions. Figure 19 depicts the geometric definitions. Mathematically. 



Figure 19. Target and Missile Flight Path Angle Definitions 
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(3.40) 


Ym.pitch ~ t 20 [^niz / ]» 

Ym_y«w~^ an [^my/^Kx]» (3*41) 

and 

Tt_piich = * an [Va/VVTVj. (3.42) 

Yi^w = tan- , [v iy /V lx ]. (3.43) 

4. Velocity and Acceleration in the Pitch Plane 

The missile is controlled in three dimensional space by generating 
acceleration commands in two orthogonal planes. These planes are defined as 
the pitch plane and the yaw plane. The yaw plane is taken as simply the 
horizontal XY plane and the pitch plane is the orthogonal plane rotated by the 
angle Oyaw. Figure 20 depicts the geometric definition of the pitch plane. Given 
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the missile velocity, and the angles a yzw and Ym_yaw. the velocity in the pitch 
plane is 

= V m COs(Ym_y«w ~ ®yaw)* (3.44) 

From the basic relationship for the missile acceleration the missile 
pitch acceleration is defined as 

i ^m_pitcb = Ym_ pilch 'Tm_pitch' (3.45) 

This acceleration vector is then broken down into Cartesian coordinate system 
components. Assuming that the pitch acceleration vector is perpendicular to 
the vertical line of sight vector between the missile and target. Figure 21 shows 
the orientation of the acceleration components. From this figure the following 
relationships are obtained 

*o_pi*di = "*(^m_pitch ‘^ pitch ) cos ^y*w» (3.46) 

ym_pildi = — (^m.pttch‘ s ^ n ^pitch (3.47) 



Figure 21. Pitch Plane Acceleration Components 





(3.48) 


^m_pitch = ^m_pitth ‘ cos ^ pitch • 

5. Velocity and acceleration in the Yaw Plane 

The missile yaw plane is depicted in Figure 22. Given the missile 
velocity and the angle Y m _pitch. the velocity in the yaw plane is simply 

YB»_y*w = ‘ C0s Ym_piich* (3.49) 



Again, from the basic missile acceleration relationship, the missile yaw 
acceleration is defined as 

31 ^n_y*w ‘ Ym_y*w* (3.50) 

The missile yaw acceleration components are depicted in Figure 23 and are 
given mathematically as 

x m_ytw = ’ShlCtytw* 


(3.51) 






y®_y*w ~ ^m_yaw' cos 


(3.52) 


yaw 



Figure 23. Yaw Plane Acceleration Components 


6. Total Missile Acceleration Components 

The input to the missile state equation is the vector consisting of the 
three Cartesian components of the overall missile acceleration vector. Given 
the missile acceleration components in both the pitch and yaw planes, the total 
missile acceleration components are 


3 x n_pitdt + x m_yaw» 

(353) 

** = Xm.pilcfa + yin_y«w* 

(354) 

and 


(355) 

Finally, the total missile acceleration is 

Am m Tl*m 2 +9m 2 +*m 2 - 

(356) 


7. Target Acceleration Components 

In this simulation the target acceleration components consist of the 
three Cartesian components of the overall target acceleration vector. These 
components are used as the input to the target kinematic equations. The input 
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vector is determined by the type of target maneuver desired. The overall target 
acceleration is simply 

A l = A /x t 2 + y t 2 + z t 2 . (3.57) 

8. Closing Velocity and Time To Go 

From equation (3.48) the velocity of the target in the pitch plane is 

simply 

^t_pitch = % COs(Yt_y*w — Oyg W ). (3.58) 

Now, the range rate, R, is found by simply projecting the missile and target 
pitch plane velocities along the vertical line of sight from the missile to the 
target Mathematically 

^ “ Vt_pitch COs(Y t pilcb — Opitch)“Vjn pitch pilch )• (3*59) 

Knowing that the range rate is the negative of the closing velocity, the time to 
go is simply 

Time - to - go = ——^— = -iL (3.60) 

^closing ~ K 

where R is the range from the missile to the target 

D. DIGITAL SIMULATION USING STATE SPACE METHOD 
1. Discrete State Equation Definition 
Given the continuous state equations 

x(t)» Ax(t)+Bu(t), (3.61) 

y(t) * Cx(t)+Du(t), (3.62) 

a more convenient form for digital simulation is 

x(k+l)*Ox(k)+ru(k), (3.63) 

y(k+l) = Cx(k)+Du(k). (3.64) 
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The discrete state matrices are defined as 

<D = e ADT = I + A • DT+(1 / 2!) A 2 • DT 2 + (1 / 3!) A 3 • DT 3 +..., (3.65) 

r = |i-DT+(1/2!)ADT 2 +(1/3!)A 2 DT 3 +...]b, (3.66) 

and for u(t) constant over [(k+l)-k]*DT. 

2. Missile Subsystem Discrete State Equations 

The missile/target scenario is simulated digitally using the MATLAB 
software package. A built in system function is used to convert the continuous 
state equations to discrete time state equations. The discrete state equations 
used in the simulation are defined as follows for the seeker head, autopilot. 


missile dynamics and target dynamics respectively; 

x sh (k+l) = O sh x sh (k)+r sh u sh (k), (3.67) 

x«p(k + l) = ^ap x ap( k ) + r ap u ap( k )* (3-68) 

x m ( k +l)*4 > m x m ( k )+r m u m (k), (3.69) 

and 

x t (k + l)-O t x t (k)+r t u t (k). (3.70) 

E. FORWARD TIME AND ADJOINT MODEL ! 


To accurately and simply assess miss distances for various target 
maneuvers, a two dimensional adjoint model is used. First, a forward time 
model is developed for the two dimensional missile/target engagement Figure 
24 shows the block diagram of the forward time model. All system inputs are 
converted to impulses for subsequent conversion to the adjoint model. 

The forward time mode is used to show the miss distance due to heading 
error (where heading error is defined as the error in missile heading from a 
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collision course) and the miss distance due to target maneuver. The missile to 
target range is defined as 

R = V c • time - to - go (3.71) 

where 

time-to-go = -t. (3.72) 

This model shows the miss distance throughout the course of the engagement. 
The last value of miss distance is the miss distance at Tnn.i. 

An adjoint model is used to analyze linear time varying systems [Ref. 6]. 
The major advantage of the adjoint technique is that the miss distance for all 
final times is generated in one run rather than the multiple runs required when 
using a Monte Carlo simulation. Following the rules for adjoint construction 
[Ref. 6], Figure 24 is converted to the adjoint model depicted in Figure 25. 
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IV. SIMULATION RESULTS 

A. OVERVIEW 

This section presents the results of the computer simulations. The scenario 
is conducted for several different initial conditions and with several different 
target maneuvers. The following assumptions are made throughout: 

1. The missile is limited to 20 g's in either plane. 

2. The acceleration due to gravity is ignored. 

3. The target is Capable of instantaneous acceleration. 

4. The target is limited to a maximum of 8.0 g's. 

5. The maximum missile speed is 3000 feet per second. 

6. The maximum target speed is 1500 feet per second. 

7. The proportional navigation constant is 4. 

8. The missile is pointed toward the target at launch. 

B. THREE DIMENSIONAL ENGAGEMENTS AGAINST 
VARIOUS TARGET MANEUVERS 

1. Scenario 1: A Constant Velocity Target 

The simulation is initially conducted with the missile engaging a 

constant velocity target Figure 26 depicts the scenario geometry. The initial 


conditions are 

xm(0) 

s, 0 feet 


yn*0) 

ss 0 feet 


zm(0) 

= 0 feet 


vmx(0) 

■ 3000 fps 
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Figure 26. Constant Velocity Target Scenario Geometry 

vmy(O) = Ofps 
vmz(0) = 0 fps 
xt(0) = 30000 feet 

yt(0) = Ofeet 

zt(0) « 500 feet 

vtx(0) = Ofjps 
vty(0) * 1000 fps 
vtz(0) = Ofps. 

Figure 27 depicts the scenario that reflects these initial conditions. 
In this case, the aircraft is attempting to stay on the outer edge of the missile 
envelope. 


38 








Figure 27. Tactical Scenario for Constant Velocity Target 
2. Scenario 2: A Constant Acceleration Target 

In this scenario, the target is accelerates at a constant rate of IS 
feet/sec/sec in the y direction. Figure 28 shows the scenario geometry. The 


initial conditions arc 

xm(b) 

= 0 feet 


ym(0) 

= 21000 feet 


zm<0) 

= Ofeet 


vmx(O) 

= 2100 fps 


vmy(0) 

= -2100 fps 


vmz(0) 

= Ofps 


xt(0) 

= 21000 feet 


yt(0) 

= Ofeet 


*(0) 

= 500 feet 
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Figure 28. Constant Acceleration Target Scenario Geometry 


vtx(O) = 0 fps 
vty(0) = 1000 fps 
vtz(0) = 0 fps. 
xddt(0) = Ofps 2 
yddtfO) * 15 fps 2 
zddt(0) = 0 fps 2 
0d(O) = -45*. 

Figure 29 depicts the scenario that reflects these initial conditions. In 
this case there is overlapping missile envelope coverage and the aircraft is 
attempting to strike a target located within this region. The missile is launched 
when the target enters the 5NM radius. 
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Figure 29. Tactical Scenario for Constant Acceleration Target 

3. Scenario 3: A Two Dimensional Target Acceleration 

Figures 26 and 27 depict the initial conditions for this scenario. The 
target performs a 6.3 g level turn toward the missile at a range of 12*000 feet. 
The target acceleration is 

x t *-6.5-32.2 sin(Y l y*w) (4.1) 

y,=6.5-32.2 cos(Y t (4.2) 

it *0.0. (4.3) 
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4. Scenario 4: A Three Dimensional Target Acceleration 

Figures 28 and 29 depict the initial conditions for this scenario. The 
target performs a three dimensional maneuver toward the missile and into the 
vertical at a range of 12,000 feet The target acceleration is 


x t =-6.5-32.2 sm(Y t 

(4.4) 

y t *6.5'32.2cos(Y tylw ) 

(4.5) 

Zj = 6.5*32.2 * cos^Yi_piich )• 

Three Dimensional Simulation Results 

(4.6) 


Figures 30-75 display the results of the three dimensional simulations. 
Figures 30-40 relate to the constant velocity target scenario. Figures 41-52 
relate to the constant acceleration scenario. Figures 53-64 relate to the two 
dimensional target acceleration scenario. Figures 65-75 relate to the three 
dimensional target acceleration scenario. 
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Figure 30. Scenario 1: Missile to Target Range 



Figure 31. Scenario 1: Missile Acceleration |A a | 
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Figure 32. Scenario 1: Missile Pitch Flight Path Angle T- 
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Figure 34. Scenario 1: Missile Yaw Flight Path Angle Ym_y*w 
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Figure 36. Scenario 1: Seeker Head Pitch Angle P piteb 



Figure 37. Scenario 1: Seeker Head Pitch Angle Rate P pi(ch 
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Figure 38. Scenario 1: Seeker Head Yaw Angle p„. 



Figure 39. Scenario 1: Seeker Head Yaw Angle Rate 
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Figure 41. Scenario 2: Missile to Target Range 
TARGET VELOCITY VS TIME 



TIME 

Figure 42. Scenario 2: Target Velocity V t 
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Figure 43. Scenario 2: Missile Acceleration |A, 



Figure 44. Scenario 2: Missile Pitch Flight Path Angle y m pitch 
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Figure 46. Scenario 2: Missile Yaw Flight Path Angle y m ^ 
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;ure 47. Scenario 2: Missile Yaw Flight Path Angle Rate Y m _yai 
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lire 57. Scenario 3: Missile Pitch Flight Path Angle Rate f«_pttch 




Figure 58. Scenario 3: Missile Yaw Flight Path Angle 7 a _j aw 
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Figure 59. Scenario 3: Missile Yaw Flight Path Angle Rate ^ 



Figure 60. Scenario 3: Seeker Head Pitch Angie p pitC || 
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Figure 63. Scenario 3: Seeker Head Yaw Angle Rate 
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Figure 64. Scenario 3: Three Dimensional Plot 
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Figure 66. Scenario 4: Missile Acceleration |A n ] 
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67. Scenario 4: Missile Pitch Flight Path Angle Y m pitch 



Figure 68. Scenario 4: Missile Pitch Flight Path Angle Rate Ymjiitcb 
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Figure 71. Scenario 4: Seeker Head Pitch Angle P pitc i, 



Figure 72. Scenario 4: Seeker Head Pitch Angle Rate P p itcb 
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Figure 73. Scenario 4: Seeker Head Yaw Angle Py** 
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Figure 75. Scenario 4: Three Dimensional Plot 
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C. MISS DISTANCE ASSESSMENT 


A two dimensional adjoint model is used to investigate the missile to target 
miss distance due to a variety of evasion maneuvers. The optimal time to 
initiate these maneuvers is determined. 

Three maneuvers are used: 

1. Step Target Acceleration. 

2. Barrel Roll. 

3. Split 'S’. 


Additional variables are 

v c 

s 

closing velocity. 

% 

= 

target acceleration magnitude. 

V m HE 

2 

missile velocity due to heading 

© t 

s 

error. 

maneuver frequency. 

L 

= 

half period of Split 'S’. 

The following initial conditions are constant throughout all simulations: 

N 

=r 

4.0. 

T f 

= ' 

6.0 seconds. 

V c 

s 

2500 feet per second. 

VmHE 

2 

0.0 feet per second. 


2 

1.0, second. 


S' 

0.1 second. 


1. Step Target Acceleration 

The step acceleration is simulated by multiplying a step input signal by 
the magnitude of the target acceleration, q t . Target acceleration.' of 1.0, 6.0, 
and 8.0 g's are simulated. 
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2. Barrel Roll 

A two dimensional Barrel Roll maneuver can be represented by a 
sinusoid of frequency G) t . A shaping filter is used to simulate the maneuver. 
[Ref. 7] Figure 76 depicts a typical Barrel Roll maneuver and Figure 77 shows 
the shaping filter equivalent where (0 t equals 1.0 radian/second. The filter is 
converted to a block diagram which is simulated using state variable design. 
The block diagram is shown in Figure 78. A 4.0 and 6.0 g Banel Roll maneuver 
are simulated. 
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« White Noise With Power 
Spectral Density 4.(0 
<^0 = H, 2 /T f OStSTf 

_ ■ 0, otherwise 

Figure 77. Shaping Filter Equivalent of a Barrel Roll 



Figure 78. Barrel Roll Block Diagram 

3. Split ’S’ 

A two dimensioned Split 'S' maneuver is represented by a periodic 
square wave with period 2L. A shaping filter is also used to simulate the 
maneuver. [Ref. 7] Figure 79 depicts a typical Split 'S' maneuver and Figure 
80 shows the shaping filter equivalent where L equals 1.0 second. The block 
diagram of this filter is shown in Figure 81. A 4.0 and 6.0 g Split 'S' maneuver 
are simulated. 












4. Forward Time and Adjoint Model Simulation Results 

Figures 82-95 show the results of the forward time and adjoint model 
simulations. Figures 82-87 show the results of the step acceleration maneuver. 
Figures 88-91 show the results of the Barrel Roll evasion maneuver. Figures 
92-95 show the results of the Split ’S’ evasion maneuver. 



Figure 79. Typical Split 'S' Maneuver 
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U 3 - White Noise With Power 
Spectral Density 4.(0 



• 0 , otherwise 


Figure 80. Shaping Filter Equivalent of Split 'S’ Maneuver 




















Figure 82. Forward time, 1 g Step Acceleration 



Figure 83. Adjoint, 1 g Step Acceleration 
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Figure 84. Forward Time, 6 g Step Acceleration 

























Figure 86. Forward Time, 8 g Step Acceleration 



Figure 87. Adjoint, 8 g Step Acceleration 
























Figure 88. Forward Time, 4 g Barrel Roll 



Figure 89. Adjoint, 4 g Barrel Roll 



























































































































V. CONCLUSIONS AND RECOMMENDATIONS 
A. CONCLUSIONS 

Clearly, the classical proportional navigation model of the missile guidance 
and control system is sufficient for simple trajectory study. The use of the 
seeker head angle rate, |3, as the estimate of the actual line of sight rate, a, 
gives accurate enough information fpr target tracking in a benign environment 
The three dimensional model of the missile tracks a maneuvering target in one, 
two, or three dimensions with enough accuracy to cause an aircraft kill. A 
navigation ratio of 4.0 generates sufficient missile acceleration commands for 
the simplified scenarios studied. 

The adjoint model proves to be a very useful indicator of miss distance 
properties. With the step acceleration evasion maneuver, the optimal time to 
initiate it is about 2.0 seconds, or 5000 feet (closing velocity is 2500 feet per 
second), prior to missile impact The miss distance increases as the target's g 
loading increases. The maximum miss distance of 140 feet is achieved with an 
8.0 g maneuver. 

With the Barrel Roll maneuver, the optimal time to initiate it is about 3.5 
seconds, or 8750 feet prior to missile impact The miss distance again 
increases with increased target g loading. The maximum miss distance of 5000 
feet is achieved with a 6.0 g maneuver. 

With the Split 'S' maneuver, the optimal time to initiate it is about 1.8 
seconds, or 4500 feet prior to missile impact. The miss distance increases with 
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increased target g loading. The maximum miss distance of 2600 feet is 
achieved with a 6.0 g maneuver. 

Oearly, the maximum miss distance and optimal time to maneuver depends 
on the evasion maneuver used. In general, it is recommended that a maximum 
g evasion maneuver be initiated when the missile is between 1.0 and 1.5 miles 
prior to impact 

B. RECOMMENDATIONS FOR FURTHER STUDY 

The simulations-developed in this thesis provide a very good generic model 
which can be easily modified for more complex missile systems. The programs 
are easily adjusted for systems with different system time constants and 
navigation ratios. 

The problems associated with measurement noise and seeker head 
measurement errors might be investigated to see their role in miss distance 
evaluation. 

It is also recommended that target Electronic Counter Measures (ECM) be 
added to the simulation of the evasion maneuvers to see their effects on miss 
distance. 
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APPENDIX A-THREE DIMENSIONAL PROGRAM 


clear 

clg 

% 

% LT Francis C. Lukenbill, USN 

% 08 August 1990 

% This program simulates a 3 dimensional target/missile 
% engagement scenario using classical proportional navigation. 

% 

% DEFINE STATES 

% 

% Missile 

% 

% ms=[xm - missile x coordinate 

% xdm - missile velocity-x direction 

% ym - missile y coordinate 

% ydm - missile velocity-y direction 

% zm - missile z coordinate 

% zdm * missile velocity-z direction] 

% 

AM=[0 10 0 0 0 
000000 
000100 
000000 
oodooi 
000000 ]; 

BM=[000 

100 

000 

010 
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000 
0 01 ]; 

% 

% Seeker Head 

% 

% beta=[beta pitch - seeker head pitch angle 

% betad pitch - seeker head pitch angle rate 

% beta yaw - seeker head yaw angle 

% betad yaw - seeker head yaw angle rate] 

% 

AS=[ 0 10 0 
-100 -20 0 0 
0 0 0 1 
0 0 - 100 - 20 ]; 

BS=[ 0 0 
100 0 
0 0 
0 100 ]; 

% 

% Autopilot 
% 

% gammadm=[gammadm pitch - pitch body angle rate 

% gammadm yaw - yaw body angle rate ] 

* ' 

AM-10 

0-1]; 

BP=fl 0 

0 1 ]; 

% 

% Target 

'% 

% t*=[xt - target x coordinate 
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% 

xdt - target velocity-x direction 

% 

yt - target y coordinate 

% 

ydt - target velocity-y direction 

% 

zt - target z coordinate 

% 

zdt - target velocity-z direction] 

% 



AT=[0 10000 
000000 
000100 
000000 
000001 
000000 ]; 

BT=[0 0 0 
100 
000 
010 
000 
00 1 ]; 

% 

% DISCRETE REPRESENTATION 

% 

dt=.01; 

[phis.de ls]= c2d(AS3S,dt); 

[phinvielm]= c2d(AM3M t dt); 
[phia3ela]= c2d(AP3P,dt); 

[phit,delt]= c2d(AT3T,dt); 

tfinal= 15.0; 
kmax=tfinal/dt +1; 

% 

% INITIALIZE VARIABLES 

% 




% navigation ratio 
% 

nr=4.0*, 

NR=[nr 0 
0 nr]; 

% 

% initial seeker head angles and angle rates 
% 

beta_pitchO =0.0; 
betad_pitch0=0.0; 
beta_yawO =0.0; 
betad_yawO =0.0; 


beta(:, 1 )=[beta_pitch0 
betad _pitchO 
beta_yaw0 
betad_yaw0 ]; 

% 

% initial missile body angle rates 
% 

gammadm_pitch0=0.0; 
gammadm_yawO =0.0; 

g:immadm(:,l)=[ganiinadm_pitchO 
gammadm_yawO ]; 


xmO: 


initial missile states 


= 0 . 0 ; 

=3000.0; 

i =0.0;; 
^|dmO=0.0; 
zm0=0.0; 
zdmOO.O; 


xdm0=3 
y|nO 
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ms(:,l)=[xmO 

xdmO 

ymO 

ydmO 

zmO 

zdmO]; 

% 

% initial target states 
% 

xtO = 30000.0; 
xdt0=0.0; 
ytO =0.0; 
ydtO= 1000.0; 
ztO =500.0; 
zdt0=0.0; 

ts(:,l)=[xtO 

xdtO 

ytO 

ydtO 

ztO 

zdtOJ; 

% 

% initial range infonnadon 
% 

rxO=(xtO-xmO); 

rx(l)=rxO; 

ryO=(ytO-ymO); 

ry(l>=ryO; 

rzO=(ztO-zmO); 

iz(l)=rzO; 

rmO=sqrt(xmO A 2 + ymO^ + zmO A 2); 
rm(l)=rmO; 





itO=sqrt(xtLi A 2 + ytO^ + ztO^); 
rt(l)=rtO; 

rO =sqrt((xtO-xmO) A 2 + (ytOymO) A 2 + (ztO-zmO)^); 
i<i)=tO; 

% 

% initial time 
% 

time(l)=0.0; 

% 

% SIMULATE THE SYSTEM 

% 

for (i=l±max-l) 

% 

% calculate missile and target speeds 

% 

vm(i)=sqrt(ms(24) A 2 + ms(4,i) A 2 + ms(6,i) A 2); 
vt(i)=sqrt(ts(24) A 2 + ts(4,i) A 2 +tsfoi)^); 

% 

% calculate line-of-sight angles 

% 

sigma_pitch(i)=atan2((ts(54)-ms(5 4)),sqrt((ts( 14)-ms( 14)) A 2. 
+{ts(34>ms(34)> A 2)); 

sigma_yaw(i)=atan2((ts(3 4)-ms(34)),(abs(ts( 14>ms( 14)));; 

sigmat:4)=isigina_pitch(i) 
sigma_yaw(i) ]; 

% 

% calculate missile and target line-of-sight angles 

% 

sigmam_pitch(i)=atan2(ms(54),sqrt(ms(14) A 2 + ms(34) A 2)); 
sigmam_yaw(i) =a:an2(ms(34),ms(14)); 
sigmat_pitch(i)=atan2<ts(54)^qn(ts( 14^2 +15(34^2)); 
*igmat_yaw(i) =atan2(ts(34),ts(14»; 
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* if* 


% 

% update seeker head states 
% 

beta(:4+l)=phis*beta(:4) + dels*sigma(: 4 >; 

% 

% set tp seeker head angle rate vector 
% 

betad(: 4 )=[beta( 24 ) 

beta(44)l; 

% 

% calculate missile and target flight path angles 
% 

gammam_pitch(i)=atan 2 (ms( 64 ).sqrt(ms( 24) A 2 + tuztf'/'-j;. 
gammam_yaw(i) =*tan 2 (ms( 44 ) 4 ns( 24 )); 
gammat_pi^h(i)=atan 2 <ts( 64 ),sqrt(ts( 24 >*2 f ts(4 i/ v 2)); 
gaminai_yaw(i) =atan2(ts(44).ts(24)); 

% 

% update autopilot states 
% 

gammadm(: 4+1 )=phii‘*gamrrwir»r 4 . •> rtrla*NR*be»ad(: 4 ); 
compute missile vekx 'vy in R*Z (pitcn,» plane (R’=Rxy) 
if (sigma_pitch(i>-= 0 . 0 : 

vmLjntch(i)^ri>s(vin(i''* , ‘:os(iigm > ‘.Vgwnmam - yaw(i))); 

% 

% limit pitch acceleration to approximately 20 g’s 

% 

if (abs(vm_pitch(i)*gammadm(l4))<s=644.0) 
am_pitch(i)=vm_pitch(i)*gammadm(l 4 ); 
else 

am_pitch(i)=644.0*sign(vm_pitcli(!)*gammadm( 1 ,i)); 
end 
else 
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vm_pitch(i)=0.0; 

am_pitch(i>=0.0; 

end 

% 

% calculate missile pitch acceleration vector components 
% 

am_pitch_xy(i)=abs(am_pitch(i)*sin(sigma_pitch(i))); 
xddm_pitch(i) =-am_pitch_xy(i)*cos(sigTna_yaw(i)); 
yddm_pitch(i) =-am_pitch__xy(i)*sin(signia_yaw(i)); 
rddm_pitch(i) =am_piU:h(i)*cos<sigma_pitch(i)); 

% 

% compute missile velocity in the XY (yaw) plane 
% 

if ((sigii-am_yaw(i) - sigmnt_yaw(i))~=0.0) 

vm_yaw(i)=abs(vm(i)*cos(gainmam_pitch(i))); 

% 

% limit yaw acceleration to approximately 20 g's 
% 

if(abs(vin_yaw(i)*gantinadni(24))<=644.0) 

am_yaw(i)=vm_yaw(i)*gammadm(24); 

else 

an\_yaw(i)=644.0*sign(vnt_yaw(i)*gaimnadm(24)); 

end 

else 

vm_yaw(i)=0,0; 

am_yaw(i>0.0; 

end 

% 

% calculate missile yaw acceleration vector components 
% 

xddm_yaw(i>=-ain_yaw(i)*sin(sigma_yaw(i)); 
yddm_yaw(i)ss am_yaw(i)*cos(sigma..yaw(i)); 

compute overall missile acceleration components 
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ip ?p ip 


% 

xddm(i)=xddm„pitch(i) + xddm_yaw(i); 
yddm(i)=yddm_pitch(i) + yddm_yaw(i); 
zddm(i)=zddm_pitch(i); 

% 

% compute total missile acceleration magnitude 
% 

am(i)=sqrt(xddm(i) A 2 + yddmCi^ + zddm(i) A 2); 

% 

% generate missile input vector 

% 

um=[xddm(i) 

yddm(i) 

zddm(i)]; 

% 

% update missile states 

% 

ms(:4+l)=y.4um*ms(:,i) + delm*um; 

% 

% set target acceleration components 

% 

if (r(i)<=0.0) 

xddt(i)=-6J*32J2*sin(gammat_yaw(i))’ 
yddt(i)= 6J*32J2*cos(gammat_yaw(i)); 
zddt(i)=0.0; 

else 

xddt(i)=0.0; 
yddt(i)=0.0; 
zddt(i)=0.0; 
end 

compute target acceleration magnitude 
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at(i)=sqrt(xddt(i)^2 + yddt(i)^2 + zddi(i)*2); 

% 

% generate target input vector 
% 

ut=[xddt(i) 

yddt(i) 

zddtCi)]; 

% 

% update target states 

% 

ts(:4+l)=phit*ts(:4) + delt*ut; 

% 

% calculate time to go 

% 

vt_pitch(i)=abs(vt(i)*cos(sigma_yaw(i)-gammat_yaw(i))); 
rdot(i)= vt_pitch(i)*cos(gamm4 k t_pitch(i)-sigma_pitch(i))... 

- vm_pitch(i)*cos(gammam_pitch(i>-sigma_pitch(i»; 
ttg(i)=T(i)/(-rdot(i)); 

% 

% calculate updated range information 

% 

r<i+l)=sqrt((ts(U+l)-ms(U+l»^2 + (tsOd+l^msOa+l)^.. 
+ (ts(54+l)-ms(54+l)) A 2); 

rm(i+l)=sqrt(ms(14+l) A 2 + ms(34+l) A 2 + ms(54+l) A 2); 

rt(i+l)=sqrt(ts(l 4+1 ^2 + ts(34+l) A 2 + ts(54+l) A 2); 

rx(i+lHts(14+l)-ms(14+l)); 

xy(i+l)=(t$(34+l) - ms(34+l)); 

rz(i+l)=(ts(54+l) - ms(54+l)); 

% 

% set up missile and target trajectory data for plotting in 3-D 
% 

missile(i,:)=[ms(14) ms(34) ms(54)l; 
target:) =[ts(14) ts(34) ts(54)l; 

% 
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% update time vector 

% 

time(i+l)=time(i) + dt; 

% 

% check to see it engagement is at closest point of approach 
% 

if(r(i)<r(i+l)),break,end 

end; 

% 

% PLOT RESULTS 

% 

% time to go 

% 

plot(r( 1 :i-1 ),ttg( 1 :i-1 )),grid,xlabelOR ange - FEET), 
ylabelCTime-to-go - SECONDS');titleCRANGE vs TIME TO GO’); 
pause,clg 
% 

% range data 

% 

plot(timej),grid^dabel(TIME , ),ylabel( , FEET) 
tiUeCRANGE VS TTME’) 
pause,clg 
% 

% missile and target velocity information 

% 

plot(time(l;i),vm),gricUlabel(TIME),ylabel(TEET/SEC) 
titleCMISSILE VELOCITY VS TIME’); 
pause,clg 

plot(dme( t ln),vt),gri(Lxlabel(TIME , ),ylabel(TEET/SEC) 
tiUe(TARGET VELOCITY VS TTME’); 
pause,clg 
% 

% missile and target acceleration information 

% 
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plot(time( 1 :i-1 ),am( 1 :i-1 )),grid,xlabel(TIME'),ylabel(TEET/SEC A 2’) 

titlcCMISSILE ACCELERATION VS TIME'); 

pause.clg 

plot(time(ln),at),grid^dabel(TIME'),ylabel('FEET/SEC A 2') 
titlcCTARGET ACCELERATION VS TIME'); 
paasc,cig 
% 

% flight path angles 
% 

plot(time(l:i),ganimain,pitch),gridptlabelCTIME'),ylabelCRADIANS’) 
titleCMSL PITCH FUGHT PATH ANGLE VS TIME'); 
pause,clg 

plot(tfjne( 1 :i),gammam_yaw),grid,xlabelCTIME'),ylabel('RADIAN S') 
titleCMSL YAW FUGHT PATH ANGLE VS TIME); 
pause,clg 

plot(tiine(l d-l),gammadm(l,(l :i- 
1 ))),grid t xlabel(TIME’),ylabel('RADIANS/SEC*) 
titleCMSL PITCH FUGHT PATH ANGLE RATE VS TIME’); 
pause,clg 

plot(time<ln-l),gainmadm(2,(ln- 
l))),grid^labelCnME'),y)abelCRADIANS/SEC) 
titleCMSL YAW FUGHT PATH ANGLE RATE VS TIME’); 
panse^lg 
% 

% seeker head angle information 
% 

plot(time( 1 :i-1 ),beta( 1,(1 :i-1 ))),gricUlabelCTIME'),ylabelCRADIANS') 

titfeCSEEKER HEAD PITCH ANGLE VS TIME'); 

panse^clg 

plot(time(l:i-l),beta(2,(l:i-l))),grid t xlabel(T[MF),ylibelCRADIANS/SEC) 

titleCSEEKER HEAD PITCH ANGLE RATE VS TIME); 

pause,clg 

plot(tune( 1 :i-1 ),beta(3,( 1 :i- l))),grid,xlabelCTIME’) f ylabelCRADIANS’) 
titf eCSEEKER HEAD YAW ANGLE VS TIME'); 
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pause, dg 

plot(time(l :i-l),beta(4,(l :i-l))),grid,xlabel(T[ME , ),ylabel( , RADL\NS/SEC) 

titleCSEEKER HEAD YAW ANGLE RATE VS TIME'); 

pause,clg 

end; 



APPENDIX B-FORWARD TIME AND ADJOINT PROGRAMS 


% 

% Lukenbill, F. C. 20 November 1990 

% This program simulates the forward time and adjoint models with a 
% step acceleration evasion maneuver. 

% 

% 

% Initialize Variables 
% 

tm=1.0; 

tshvO.l; 

a=l/tm; 

b=2*(l/tsh); 

css(lAsh)*(l/tsh); 

N=4.0; 

Vc=2500; 

Tfc=6.0; 

ki=l/(Vc*Tf); 

HE=0.0*pi/180; 

Vm=2000.0; 

VmHE=Vm*sin(HE); 

Atgt=257.6;; 

% 

% Input State Matricies 
% 

A=[0 1 0 0 0 0 0 0 -VmHE 

0 0 1 0 0 0 0 0 0 

0 0 -a 0 a*N*Vc 0 0 0 0 

0 0 0 0 1 0 0 0 0 

-c*lri 0 0 -c -b c*ki 0 0 0 

0 0 0 0 0 0 1 0 0 
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0 0 0 0 0 0 0 Atgt 0 

0 00 0 00000 
0 0 0 0 0 0 0 0 0 ]; 

% 

B=[0 

0 

0 

0 

0 

0 

0 

1 

l]; 

% 

C=[-l 0 0 0 0 1 0 0 0]; 

% 

% Input Timing Information 
% 

dt=0.01; 

kmax-Tf/dt+1; 

% 

% Set up impulse function 
% 

impulse=zeros( 1 Jonax); 
impulse(l)=l/dt; 

% 

% Initialize vectors 
% 

x=zerosf Janax); 
y=zeros(lJanax); 
time=zeros( 1 ,kmax); 

% 

% Simulate the System (Forward Time Solution) 

% 
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for (i=l:kmax-l) 

ki=l/(Vc* v Tf-dme(i)))r 

A=[0 1 0 0 0 0 0 0 -VmHE 

0 01 0 0 0 0 0 0 
0 0 -a 0 a*N*Vc 0 0 0 0 

0 0 0 010 0 0 0 
-c*ki 0 0 -c -b c*ki 0 0 0 

0 0 0 0 0 0 1 0 0 

0 0 0 0 0 0 0 Atgt 0 

0 0000000 0 
0 0 0 0 0 0 0 0 0 ]; 

[Phifwd J)elfwd] =c2d( A,B ,dt); 
x(:4+l) * Phifwd*x(:4) + Detfv'd*impulse(i); 
y(:4+l) = C*x(:4+1); 
tiine(i+l )=time(i)+dt; 
end; 

y(i) 

% 

% Plot Results 
% 

plot(time(l:i),y(l,(l:i»);title< , Fofward Tune - 8 g Step Acceleration’); 
xlabelOIme’fc.ylabelOy miss’);grid;pause; 

% 

% Reinitialize variables and vectors 
% 

dt=0.01; 

kmax-Tf/dt+1; 

impulse=zeros( 1 Jcmax); 

impulse(l)=l/dt; 

x=zeros(9,kmax); 

y=zeros(ljonax); 

time=zeros(14anax); 

time(l)=dt; 

% 
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% Simulate the System (Adjoint Solution) 

% 

for (i=l:kmax-l) 

ki=l/(Vc*(tir.ie(i))); 

A-[0 1 0 0 0 0 0 0 -VmHE 

0 0 1 0 0 0 0 0 0 

0 0 -a 0 a*N*Vc 0 0 0 0 

0 0 0 0 1 0 0 0 0 

-c*ld 0 0 -c -b c*ki 0 0 0 

0 0000C10 0 

0 0 0 0 0 0 0 Atgt 0 

0 0 0 0 0 0 0 0 0 

0 0 0 0 0 0 0 0 0 ]; 

[Phiadj,Deladjl=c2d(A',C,dt); 
x(:,i+l) = Phiadj*x(:,i) + Deladj*impu!se(i); 
y(:4+l) *= B'*x(:,i+1); 
time(i+l )=time(i)+dt, 
end; 

><i) 

% 

% Plot Results 
% 

plotCtuneClrO.yCl.vla^ititleCAdjoint * 8 g Step Acceleration'); 
xlabe!(Time');,ylabel(’y miss’);giid;pause; 



% 

% Lukenbill, F. C. 20 November 1990 
% This program simulates the forward time and adjoint models with a 
% Barrel Roll maneuver. 


% 

% 

% Initialize Variables 


O', 


tm-1.0; 

tsh=0.1; 

a=l/tm; 

b=2*(l/tsh); 

c=(l/tsh)*(l/tsh); 

N=4.0; 

Vc=2500; 

Tf=6; 

ki=VfVc*Tf); 

HE=u.0*pi/180; 

Vm=2000.0; 

VmHE=Vm*sin(HE); 

Atgt=193.2; 

wt=1.0; 

Vr =wt*wt; 

ATGT=wt*(Atgt*Atgt)/Tf; 

K 

% Input State Matricies 
% 

A=[ 0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 -a 0 a*N*Vc 0 

0 0 0 0 1 0 

-c*ld 0 0 -c -t c # ki 

0 G O 0 0 0 

0 0 0 0 0 0 


0 

0 

0 

0 

0 

1 

0 


0 -VmHE 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 0 0 
0 0 1 0 
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0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

% 

B=[0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

01 ; 

% 

C=[-l 0 0 0 0 1 0 0 0 0 0]; 
% 

% Input Timing Information 
% 

dt=0.01; 
kmax=Tf/dt+1; 

% 

% Set tip impulse function 
% 

impulse=zeros( 1 Jonax); 
impulse(l)=l/dt; 

% 

% Initialize vectors 
% 

x=zeros(l ljemax); 
y=zeros(l r kmax); 
time=zeros( 1 Janax); 


0 

0 

0 

0 


0 0 0 

0 0 0 

0 0 0 

0 ATGT 0 


0 0 
C 0 
0 1 
-W 0] 
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% 

% Simulate the System (Forward Time Solution) 
% 


for (i=l:kmax-l) 
lri= 

A=[ 0 1 

=l/(Vc*(Tf-time(i»); 

0 0 0 

0 

0 

0 

-VmHE 0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

•a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

-c*lri 

0 

0 

-c 

-b 

c*ld 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

C 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

ATGT 0 

-V 

[PhifvrdJDelfwd]=c2d(A t B,dt); 

x(:,i+l) = Phifwd*x(:4) + Delfwd*impulse(i); 





y(:,i+l) - C*x(:j+1); 
time(i+l )=time(i)+dt; 
end; 

ytf) 

% 

°h Plot Results 
% 

plot(time(l :i),y( 1X1 n)));title( , Forward Time - 6 g Band Roll'); 
xlabel(Time');,ylabelCy miss’);grid;pause; 

% 

% Reinitialize variables and vectors 
% 

dt=0.01; 
kmax=Tf/dt+1; 
impulse=zeros(lJcmax); 
impulse(l)=l/dt; 







x=zeros(lljonax); 
y=zcros( 1 Jonax); 
timc=zcros< 1 Jonax); 
time(l)=dt; 

% 

% Simulate the System (Adjoint Solution) 
% 


for (i=l Jonax-1) 
ld= 

:l/(Vc*(time(i))); 







A=[ 0 

1 

0 

0 

0 

0 

0 

0 - 

VmHE 0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

-c*ld 

0 

0 

-c 

-b 

c*ld 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 0 0 0 0 
[Phiadj,Deladj]=c2d(A\C\d:); 

0 

0 

ATGT 

0 

-w 

01; 


x(:J+l) = Phiadj*x(:J) + Deladj*impulse(i); 
y(: j+1)» B'*x(:J+l); 
time(i+l )=time(i)4dt; 
end; 

*0 

% 

% Plot Results 
% 

plot(time(l:i),y(l,(ln)));titleCAdjoint - 6 g Barrel Roll'); 

xlabelCTime’);,ylabelCy miss');grid;pause; 

end; 
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% 

% Lukenbill, F. C. 21 November 1990 

% This program simulates the forward time and adjoint models with a 
% Spilt 'S’ maneuver 
% 

% Initialize Variables 
% 

tm=1.0; 

tsh=0.1; 

a=tm; 

b=2*(l/tsh); 

c=(lAsh)*(l/tsh); 

N=4.0; 

Vc=2500; 

Tf=£; 

ld=l/(Vc*Tf); 

HE=0.0*pi/180; 

Vm=2000.0; 

VmHE=Vm*sin(KE); 

Atgt=193.2; 

L=1.0; 

al=pi A 2/L A 2; 

i2=pi*l 6* Atgt A 2/(L*pi A 2*Tt); 

«3*9*al; 

a4=a2; 

% 

% Initialize State Matricies 


% 


A=[ 0 

1 

0 

0 

0 

0 

0 

0 

-VmHE 0 

000 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

•a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

000 

-c*ld 

G 

0 

•c 

-b 

c*ld 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

GOO 
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0 000000001 010 

0 0 0 0 0 0 0 0 0 0 000 

0 000000000 000 

0 0 0 0 0 0 0 0 0 0 100 

0 0 0 0 0 0 0 ?2 0 -al 000 

0 000 000000 001 
0 0 0 0 0 0 0 a4 0 0 0-a3 0] 

% 

B=[0 

0 

0 

0 

0 

0 

0 

1 

1 

0 

0 

0 

01 ; 

% 

C=[-l 000010 0 0 0 00 0J; 

% 

% Input Timing Information 
% 

dt=0.01; 

kmax=Tf/dt+l; 

% 

% Set up impulse function 
% 

impulse=zeros( ljcmax); 
impulsc(l)=l/dt; 

% 
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% Initialize vectors 
% 

x=zeros( 13,kmax); 
y=zeros(l,kmax); 
time=zeros( 1 Jonax); 

% 

% Simulate the System (Forward Time Solution) 
% 


for (i=l±max-l) 



ld= 

:l/(Vc*(Tf-time(i))); 







> 

II 

r—* 

o 

1 

0 

0 

0 

0 

0 

0 

-VmHE 

0 

000 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

•a 

0 

a*N*Vc 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

000 

-c*ld 

0 

0 

-c 

-b 

c*ld 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

010 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

100 

0 

0 

0 

0 

0 

0 

0 

a2 

0 

-al 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

001 

0 

0 

0 

0 

0 

0 

0 

a4 

0 

0 

0 -a3 0]; 


[Phifwd,Delfwd]=c2d(A,B,dt); 
x(:4+l) * Phifwd*x(:4) + De!fwd*impulse(i); 
y(:,i+l) = C*x(:4+l); 
time(i+l)=time(i)+dt; 
end; 

>0) 

% 

% Plot Results 
% 

plot(time( 1 :i),y( 1,(1 n)));title(Torward Time 6 g Split S'); 
xlabel(Tlme');,ylabel('y miss');grid;pause; 






% 

% Reinitialize variables and vectors 
% 

dt=0.01; 
kmax=Tf/dt+l; 
impulse=zeros( 1 Jcmax); 
impulse(l)=l/dt; 
x=zeros(13Jcmax); 
y=zeros(l,kmax); 
dme=zeros( 1 Jon ax); 
time(l)=dt; 

% 

% Simulate the System (Adjoint Solution) 
% 

for (i=ldanax-l) 


ki=l/(Vc*(time(i))); 


> 

II 

f—» 

o 

1 

0 

0 

0 

0 

0 

0 

-VmHE 0 

000 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

-a 

Q 

a*N*Vc 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

000 

-c*ld 

0 

0 

-c 

-b 

c*ld 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

,1 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

010 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

100 

0 

0 

0 

0 

0 

0 

0 

a2 

0 

•al 

000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

001 

0 

0 

0 

0 

0 

0 

0 

a4 

0 

0 

0 -a3 0] 


[Phiadj,Deladj]=c2d(A\C,dt); 
x(:4+l) = Phiadj*x(;,i) + Deladj*impulse(i); 
yC4+l) = B’*x(:,i+1); 
time(i+l)=tiine(i)+dt; , 
end; 




% 

% Plot Results 
% 

plot(time( 1 d),y( 1,(1 :i)));titleCAdjoint - 6 g Split S’); 
xlabel(Time');,ylabel(’y miss');grid;pause; 
end; 





APPENDIX C-DISSPLA PLOTTING PROGRAM 
CC Lukenbill, F. C. 6/15/90 

CC This code is used to generate the three dimensional plots 
CC 

DIMENSION MSLX(300), MSLY(300),MSL^300), 
TGTX(300),TGTY(300),TGTZ(300) 

NUM=100 
DO 10 I=1»NUM 

READ (10,1000) MSLX(I),MSLY(I),MSLZ(I) 

READ (11,1000) TGTX(I),TGTY(I),TGTZ(I) 

10 CONTINUE 

1000 FORMAT (3(1X,IE15.7)) 

CALLCOMPRS 
CALL AREA2D(8.,8.5) 

CALL X3NAMECX AXIS$\100) 

CALL Y3NAMECYAXIS$\ 100) 

CALL Z3NAMECZ AXIS$M00) 

CALL VOLM3D(7.,7.,7.) 

CALL VIEW(80000.,24000.,2000.) 

CALL GRAF3D(0.,4000.,24000.,0.,4000.,24000 M 0.,400.,2400.) 
CALLGRFrn(0.,0.,0.,0.,l.,0.,0.,l.,l.) 

CALL AREA2D(7.,7.) 

CALLGRAF(0.,1.,6.,0.,1.,6.) 

CALL GRffKW) 

CALL END3GR(0) 

CALL GRFrn(0.,0.,0.,0.,l.,0.,l.,l.,0.) 

CALL AREA2C v /'.,7.) 

CALL GRAF(0.,4000,^4000,0.,4000.,24000.) 

CALL DASH 

CALL CURVE (MSLY,MSLX,NUM,0) 










CALL RESETODASH) 

CALL CURVE (TGTY,TGTX,NUM,0) 

CALL GRED(2,2) 

CALL END3GR(0) 

CALL DASH 

CALL CURV3D(MSLX>1SLY^SLX,NUM,0) 
CALL RESETCDASH*) 

CALL CURV3D(TGTX,TGTY,TGTX,NUM,0) 
CALL ENDPL(O) 

CALL DONEPL 

STOP 

END 
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